# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde
from matplotlib.patches import Ellipse
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from pyproj import Transformer
from sklearn.cluster import DBSCAN
from libpysal.weights import DistanceBand
from esda.moran import Moran
import geopandas as gpd
from shapely.geometry import Point

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False


# ==================== 1. 数据加载与预处理 ====================
def load_data():
    """加载并清洗西安POI数据"""
    try:
        df = pd.read_csv('warehouse_xian_poi.csv')
        print(f"原始数据量: {len(df)}条")

        # 数据清洗
        df = df.dropna(subset=['lng', 'lat', 'type'])
        print(f"去除空值后: {len(df)}条")

        # 筛选旅游相关类别
        tourism_types = ['scope']  # 根据实际数据调整
        df_tourism = df[df['type'].isin(tourism_types)].copy()
        print(f"筛选旅游类别后: {len(df_tourism)}条")

        # 验证坐标范围（西安主城区）
        df_tourism = df_tourism[
            (df_tourism['lng'] > 108.7) & (df_tourism['lng'] < 109.3) &
            (df_tourism['lat'] > 33.9) & (df_tourism['lat'] < 34.5)
            ]
        print(f"最终有效数据量: {len(df_tourism)}条")

        # 转换为GeoDataFrame
        geometry = [Point(xy) for xy in zip(df_tourism['lng'], df_tourism['lat'])]
        gdf = gpd.GeoDataFrame(df_tourism, geometry=geometry, crs="EPSG:4326")

        return gdf
    except Exception as e:
        print(f"数据加载出错: {str(e)}")
        return gpd.GeoDataFrame()


# ==================== 2. 空间分析核心功能 ====================
def spatial_analysis(gdf):
    """执行完整的空间分析流程"""
    if gdf.empty:
        return None, None, None

    # 坐标转换（WGS84转UTM Zone 49N）
    transformer = Transformer.from_crs(4326, 32649)
    coords = np.array([transformer.transform(point.y, point.x)
                       for point in gdf.geometry])
    gdf['x'], gdf['y'] = coords[:, 0], coords[:, 1]

    # 1. 标准差椭圆计算
    def calc_ellipse(x, y):
        cov = np.cov(x - x.mean(), y - y.mean())
        eigvals, eigvecs = np.linalg.eig(cov)
        major = 2 * np.sqrt(eigvals[0])
        minor = 2 * np.sqrt(eigvals[1])
        angle = np.degrees(np.arctan2(*eigvecs[:, 0][::-1]))
        return major, minor, angle

    major, minor, angle = calc_ellipse(gdf['x'], gdf['y'])
    center = gdf.geometry.unary_union.centroid
    center_lng, center_lat = center.x, center.y

    # 2. 核密度估计
    kde = gaussian_kde(np.vstack([gdf.geometry.x, gdf.geometry.y]))
    gdf['density'] = kde(np.vstack([gdf.geometry.x, gdf.geometry.y]))

    # 3. DBSCAN空间聚类
    coords_rad = np.radians(np.column_stack([gdf.geometry.y, gdf.geometry.x]))
    dbscan = DBSCAN(eps=0.002, min_samples=5, metric='haversine').fit(coords_rad)
    gdf['cluster'] = dbscan.labels_

    # 4. 空间自相关分析（使用投影坐标）
    try:
        w = DistanceBand.from_dataframe(gdf.to_crs(epsg=32649), threshold=2000)  # 2km阈值
        moran = Moran(gdf['density'].values, w)
        print(f"Moran's I: {moran.I:.3f} (p-value: {moran.p_sim:.3f})")
    except Exception as e:
        print(f"空间自相关分析失败: {str(e)}")
        moran = None

    return gdf, (major, minor, angle, center_lng, center_lat), moran


# ==================== 3. 增强可视化 ====================
def plot_analysis_results(gdf, ellipse_params, moran=None):
    """绘制综合空间分析结果"""
    if gdf.empty:
        return

    # 创建地图画布
    fig = plt.figure(figsize=(18, 12))
    proj = ccrs.PlateCarree()

    # 1. 主图：空间分布与椭圆
    ax1 = fig.add_subplot(221, projection=proj)
    ax1.set_extent([108.85, 109.15, 34.15, 34.45], crs=proj)

    # 地图底图
    ax1.add_feature(cfeature.LAND.with_scale('50m'), facecolor='#f5f5f5')
    ax1.add_feature(cfeature.RIVERS.with_scale('50m'), edgecolor='#66b3ff')
    ax1.gridlines(draw_labels=True, linestyle='--', alpha=0.5)

    # 绘制核密度热力图
    xi, yi = np.mgrid[108.85:109.15:200j, 34.15:34.45:200j]
    zi = gaussian_kde(np.vstack([gdf.geometry.x, gdf.geometry.y]))(np.vstack([xi.flatten(), yi.flatten()]))
    ax1.contourf(xi, yi, zi.reshape(xi.shape), levels=15, cmap='RdYlBu_r', alpha=0.6, transform=proj)

    # 绘制标准差椭圆
    major, minor, angle, clng, clat = ellipse_params
    ellipse = Ellipse(
        (clng, clat), width=minor / 50000, height=major / 50000, angle=angle,
        edgecolor='#e63946', facecolor='none', linewidth=2, transform=proj, label='Std Ellipse'
    )
    ax1.add_patch(ellipse)
    ax1.plot(clng, clat, 'ro', markersize=8, transform=proj)

    # 2. 子图1：聚类结果
    ax2 = fig.add_subplot(222, projection=proj)
    ax2.set_extent([108.85, 109.15, 34.15, 34.45], crs=proj)
    ax2.set_title('DBSCAN Clustering Results', pad=15)

    # 按聚类分组绘制
    for cluster in gdf['cluster'].unique():
        if cluster == -1:  # 噪声点
            color = 'gray'
            alpha = 0.3
        else:
            color = plt.cm.tab20(cluster % 20)
            alpha = 0.7
        cluster_data = gdf[gdf['cluster'] == cluster]
        ax2.scatter(cluster_data.geometry.x, cluster_data.geometry.y, c=[color],
                    s=15, alpha=alpha, transform=proj, label=f'Cluster {cluster}')

    # 3. 子图3：参数表格
    ax3 = fig.add_subplot(223)
    ax3.axis('off')
    params = [
        ["Center", f"{clng:.3f}°E, {clat:.3f}°N"],
        ["Major Axis", f"{major / 1000:.1f} km"],
        ["Minor Axis", f"{minor / 1000:.1f} km"],
        ["Orientation", f"{angle:.1f}°"],
        ["Area", f"{np.pi * major * minor / 1e6:.1f} km²"],
        ["Clusters", f"{len(gdf['cluster'].unique())} (incl. noise)"]
    ]

    if moran:
        params.append(["Moran's I", f"{moran.I:.3f} (p={moran.p_sim:.3f})"])

    table = ax3.table(cellText=params, loc='center', cellLoc='left', colWidths=[0.3, 0.7])
    table.auto_set_font_size(False)
    table.set_fontsize(12)
    table.scale(1, 2)

    # 全局设置
    plt.suptitle('Xi\'an Tourism POI Spatial Analysis', fontsize=16, y=0.95)
    plt.tight_layout()
    plt.savefig('xian_poi_full_analysis.png', dpi=300, bbox_inches='tight')
    plt.show()


# ==================== 主执行流程 ====================
if __name__ == "__main__":
    # 1. 数据加载
    gdf = load_data()

    if not gdf.empty:
        # 2. 空间分析
        gdf, ellipse_params, moran = spatial_analysis(gdf)

        # 3. 可视化
        plot_analysis_results(gdf, ellipse_params, moran)

        # 保存分析结果
        gdf.to_file('xian_poi_analysis_results.geojson', driver='GeoJSON')
    else:
        print("执行失败：无有效数据")